1

ECEN 4532: DSP Laboratory

Professor Michael Perkins

Lab 5

March 5th, 2018

2

Contents

1. Introduction 3

2. The JPEG Compression Algorithm 3

3. Conclusion 27

A. MATLAB Code 27

A.1. MATLAB code used to compile report 27

A.2. MATLAB functions 29

A.2.1 dctmgr.m 29

A.2.2 idctmgr.m 30

A.2.3 quantization.m 30

A.2.4 inverseQuantization.m 30

A.2.5 zigZag.m 31

A.2.6 inverseZigZag.m 31

A.2.7 run_bits_value.m 33

A.2.8 inverse_run_bits_value.m 34

3

1 Introduction

In this lab we are implementing the JPEG algorithm used for image compression. We do this so

that images can be transferred and stored efficiently and to save space. Uncompressed images

are very difficult to store, so storing these images digitally is very advantageous.

We are using an encoding/decoding algorithm in order to achieve compression on luminance

(gray-level) images. The main part of this lab is using the Discrete Cosine Transform (DCT) to

describe an image using a smaller set of coefficients to represent the original image. These values

are put through a quantization table to decrease the information needed to represent the

coefficients created through the DCT. Then we decode the information from these steps to get a

new image that is compressed that uses less space than the original picture.

2 The JPEG Compression Algorithm

To implement the JPEG compression algorithm we start by taking the forward and inverse DCT. To

implement the forward DCT the original image is divided into non overlapping blocks size 8x8.

Each block gets transformed using the 2D DCT expressed as,

These coefficients are an approximation and do introduce distortion and slight loss in image quality.

Given the coefficients 𝑓

𝑢,𝑣

̃

one can use the inverse DCT to reconstruct an estimate of the block using

4

1) The function dctmgr is included at the end of the report in section A.2.1.

2) The function idctmgr is included at the end of the report in section A.2.2.

The next step in the JPEG standard is quantization. The effect of the DCT is to create many small

coefficients that are close to zero, and a small number of large coefficients. The goal of

5

quantization is to represent a continuum of values with a finite small set of symbols. In JPEG the

quantization of the DCT coefficients other than the DC coefficient is performed as follows.

After using the quantization, we take the inverse that is defined by the map

We notice that by doing this, the compression becomes lossy during the quantization step as it

creates many zero coefficients as the loss factor becomes larger. As a result, it becomes possible

to spend fewer bits.

3) The functions quantization.m and inverseQuantization.m are using to implement this part of the

lab and are included at the end of the report in section A.2.3 and A.2.4.

4) Below are the six test images reconstructed with a loss factor of 1, 10, and 20 as well as plots of

the reconstruction error, where areas that are more yellow signify a larger error.

6

Figure 1: reconstruction and reconstruction error of barbara.pgm at LF 1

7

Figure 2: reconstruction and reconstruction error of clown.pgm at LF 1

8

Figure 3: reconstruction and reconstruction error of fingerprint.pgm

9

Figure 4: reconstruction and reconstruction error of mandril.pgm with LF 1

10

Figure 5: reconstruction and reconstruction error of roof.pgm with LF1

11

Figure 6: reconstruction and reconstruction error of straw.pgm with LF 1

12

Figure 7: reconstruction and reconstruction error of barbara.pgm with LF 10

13

Figure 8: reconstruction and reconstruction error of clown.pgm with LF 10

14

Figure 9: reconstruction and reconstruction error of fingerprint.pgm with LF 10

15

Figure 10: reconstruction and reconstruction error of mandril.pgm with LF 10

16

Figure 11: reconstruction and reconstruction error of roof.pgm with LF 10

17

Figure 12:reconstruction and reconstruction error of straw.pgm with LF 10

18

Figure 13: reconstruction and reconstruction error of barbara.pgm with LF 20

19

Figure 14: reconstruction and reconstruction error of clown.pgm with LF 20

20

Figure 15: reconstruction and reconstruction error of fingerprint.pgm with LF 20

21

Figure 16: reconstruction and reconstruction error of mandril.pgm with LF 20

22

Figure 17: reconstruction and reconstruction error of roof.pgm with LF 20

23

Figure 18: reconstruction and reconstruction error of straw.pgm with LF 20

24

In figures 1, 7, and 13 the reconstruction is very poor especially around the person as the loss

factor increases. The reconstruction error is quite large for a loss factor of 20, but I cannot really tell the

difference between loss factor of 1 and 10 as I can’t really tell different shades of blue very well, but I

decided not to scale these images to tell the reconstruction error better as I don’t see the value in it as

the difference between the images doesn’t seem to large between loss factors of 1 and 10.

In figures 2, 8, and 14 the story is similar to barbara.pgm as the image has pixels next to each

other are pretty different so the reconstruction is pretty bad for loss factor of 20. For loss factors of 1

and 10 the story is similar to babara.pgm.

In figures 3, 9, and 15 the reconstruction is probably the worst out of all the images as the loss

factor goes up. It is very clear that there are many patches of the fingerprint that are reconstructed

poorly which makes sense as fingerprints have very unique patterns.

In figures 4, 10, and 16 the reconstruction is good around the center of the mandril’s face, but

the features around the face are reconstructed poorly as loss factor goes up.

In figures 5, 11, and 17 the reconstruction on the top of the roofs is poor as loss factor goes up.

This makes sense as there are many ridges in the roof that causes the reconstruction to be horrible.

In figures 6, 12, and 18 the reconstruction seems to work reliably and there is only poor

reconstruction for small amounts of pixels rather than a large patch of pixels compared to the other

images when loss factor was high.

5) Below in figure 19 is the PSNR vs. Loss Factor for each of the 6 images. The plots have

logarithmic decay as loss factor increases which makes sense as the reconstruction gets worse

as the loss factor increases and the PSNR has a factor of log 10 in its calculation. It also makes

sense as there is significantly more noise as the reconstruction gets worse which causes the

PSNR to decrease which can be seen in the previous graphs for part 4 as loss factor increases.

25

Figure 19: PSNR vs. Loss Factor

The next part of the lab involves computing the maximum and minimum values of any DCT coefficient or

DC prediction residual can assume. When a quantized coefficient is a zero, the coefficient is not coded.

Instead, we count the number of zeros separating the non-zero coefficient from the previous non-zero

coefficient, we precede it with a code word that encodes the number of zeros between this coefficient

and the previous non-zero coefficients. To achieve this we use arithmetic encoding, and to create the

input to the arithmetic encoder we process the image’s coefficients into a triplet of numbers: nZeroes,

nBits, and value. These values are based on the DC coefficients of the image and are implemented in the

next assignment.

26

6a) Implementation of the pre-processor run_bits_value can be found at the end of the report in section

A.2.7.

6b) Implementation of the post-processor inverse_run_bits_value is at the end of the report in section

A.2.8.

After implementing the pre-processor, we go to the arithmetic encoding step by feeding in the symb

vector to the arithmetic encoder. This is done in the next assignment.

27

7) The symb vector was produced from run_bits_value and I put it into the irun_bits_value function, but

I did not have time to evaluate the performance of the arithmetic encoder/decoder.

8) This part was skipped as I ran out of time to evaluate the performance of my functions.

3 Conclusion

In this lab, by using the DCT and arithmetic encoding and decoding, we were able to implement

the JPEG algorithm. It could be seen that with greater loss factors that the image reconstruction

becomes worse. By comparing loss factor to the PSNR of the reconstructed images to the original

images as well as the compression ratio to the PSNR, we were able to visually see how well

reconstruction was for certain loss factors.

A MATLAB Code

A.1 Code used to compile report

%This is the main script for lab 5

close all;

clear all;

clc;

pathName = 'C:\Users\house_000\Documents\MATLAB\lab5\images';

fileList = dir(fullfile(pathName,'*.pgm'));

fileName = dir('*.pgm');

numberOfImages = size(fileList,1);

28

for i = 1:numberOfImages

%questions 1-4 parts:

f = imread(fileName(i).name);

figure

imagesc(f, [0 255])

colormap gray

axis square

axis off;

lossFactor = 10;

g = im2double(f);

[coefficientMatrix,temp] = dctmgr(g,lossFactor);

outputImage = idctmgr(coefficientMatrix,temp,lossFactor);

outputImage = im2uint8(outputImage);

figure

imagesc(outputImage, [0 255])

colormap gray

axis square

axis off;

reconstructionError = f - outputImage;

figure

imagesc(outputImage, [0 255])

colormap gray

axis square

axis off;

title(strcat('Reconstructed Image with Loss Factor of: ','

',int2str(lossFactor)));

saveas(gcf,strcat('Image',int2str(i),'

','lossFactor',int2str(lossFactor),'.png'));

figure,imagesc(reconstructionError);

title(strcat('Reconstruction error with Loss Factor of ', '

',int2str(lossFactor)));

saveas(gcf,strcat('Reconstruction error',int2str(i),'

','lossFactor',int2str(lossFactor),'.png'));

end

%question 5

PSNRMatrix = zeros(50,6);

parfor i = 1:numberOfImages

f = imread(fileName(i).name);

%f = imread('mandril.pgm');

g = im2double(f);

for lossFactor = 1:50

lossFactor = 5;

[coefficientMatrix,temp] = dctmgr(g,lossFactor);

symb = run_bits_value(coefficientMatrix);

%minval = min( symb );

%symb = symb - (minval - 1);

%maxval = max( symb );

%histvals = histcounts(symb, maxval);

%histvals = histvals + 1;

%bitstream = arithenco(symb, histvals);

%total_bits = length(bitstream);

%symb_length = length(symb);

%recon_symb = arithdeco(bitstream, histvals, symb_length);

29

%recon_symb = cast(recon_symb, 'double') + double(minval -1);

%irbv_coeffs = irun_bits_value(recon_symb);

outputImage = idctmgr(coefficientMatrix,temp,lossFactor);

outputImage = im2uint8(outputImage);

reconstructionError = f - outputImage;

[peaksnr,snr] = psnr(outputImage, f);

PSNRMatrix(lossFactor,i) = peaksnr;

end

figure

lossFactor = 1:50;

fileName = {'barbara', 'clown', 'fingerprint', 'mandril', 'roof', 'straw'}

for i = 1:6

subplot(2,3,i)

plot(lossFactor,PSNRMatrix(:,i));

title(strcat(fileName(i), ' ', 'PSNR vs. Loss Factor'));

end

saveas(gcf, 'LossFactorvsPSNR.png');

A.2 MATLAB functions

A.2.1. dctmgr.m

function [coefficients, temp] = dctmgr(image, lossFactor)

%dctmgr Creates JPEG coefficients for an image

%Takes a luminance (gray-level) image as an input, divives it into

%non overlapping 8x8 blocks, and uses the DCT transform on each block

according

%F(u,v) = 1/4*C_u*C_v sum(f(x,y)*cos(((2x+1)*u*pi)/16) *

%cos(((2y+1)*v*pi)/16), u,v = 0,...,7. For a given block b,

%the coefficients of the column coefficients(:,b) are organized in a

%zig-zag pattern. The zero-frequency coefficients, F(1,1), for each block

%is encoded using differential encoding: coefficients(1,1) = F_1(1,1);

%coefficients(1,2) = F_2(1,1) - F_1(1,1); where F_i(1,1) is the zero

%frequency DCT coefficient of block i

[x,y] = size(image);

index = 1;

coefficients = zeros(64,x*8);

dctMatrix = zeros(x,y);

for x2 = 1:8:x

for y2 = 1:8:y

temp = image(x2:x2+7,y2:y2+7);

firstDct = dct2(temp);

firstDct = quantization(firstDct,lossFactor);

dctMatrix(x2:x2+7,y2:y2+7) = firstDct;

%coefficient matrix:

coefficients(:,index) = zigZag(firstDct);

index = index + 1;

end

temp = coefficients(1,:);

for i = 2:y*8

coefficients(1,i) = temp(1,i) - temp(1,i-1);

end

30

A.2.2 idctmgr.m

function [output] = idctmgr(coefficientMatrix,temp,lossFactor)

%(a) It takes a matrix coeff of size 64 x Nblocks, where Nblocks is the

number of 8 x 8 blocks

% inside the image and reconstruct a luminance image.

% (b) For each given block b, the coefficients in the column coeff(:,b) are

used to reconstruct

% the block according to (3).

firstiDct = zeros(8,8);

output = zeros(512,512);

[x,y] = size(coefficientMatrix);

coefficientMatrix(1,:) = temp;

index = 1;

for x2 = 1:8:505

for y2 = 1:8:505

firstiDct = inverseZigZag(coefficientMatrix(:,index));

firstiDct = inverseQuantization(firstiDct,lossFactor);

output(x2:x2+7,y2:y2+7) = idct2(firstiDct);

index = index + 1;

end

A.2.3 quantization.m

function [output] = quantization(f_uv,lossFactor)

%Implement the quantizer

output = zeros(8,8);

table_Q = [16 11 10 16 24 40 51 61;

12 12 14 19 26 58 60 55;

14 13 16 24 40 57 69 56;

14 17 22 29 51 87 80 62;

18 22 37 56 68 109 103 77;

24 35 55 64 81 104 113 92;

49 64 78 87 103 121 120 101;

72 92 95 98 112 100 103 99];

table_Q = table_Q/255;

for u = 1:8

for v = 1:8

output(u,v) = round((f_uv(u,v)./(lossFactor.*table_Q(u,v))));

end

A.2.4 inverseQuantization.m

function [output] = inverseQuantization(f_uv,lossFactor)

%Inverse quantization:

output = zeros(8,8);

31

table_Q = [16 11 10 16 24 40 51 61;

12 12 14 19 26 58 60 55;

14 13 16 24 40 57 69 56;

14 17 22 29 51 87 80 62;

18 22 37 56 68 109 103 77;

24 35 55 64 81 104 113 92;

49 64 78 87 103 121 120 101;

72 92 95 98 112 100 103 99];

table_Q = table_Q/255;

for u = 1:8

for v = 1:8

output(u,v) = f_uv(u,v).*lossFactor.*table_Q(u,v);

end

A.2.5 zigZag.m

function [ outputMatrix ] = zigZag( inputMatrix )

index = reshape(1:numel(inputMatrix),size(inputMatrix));

%flip the matrix left-right. Then use spdiags

%spdiags will get extract the diagional elements and then

%flip the element left-right

index = fliplr(spdiags(fliplr(index)));

%reverse order of odd columns and flip up-down

index(:,1:2:end) = flipud(index(:,1:2:end));

index(index==0) = [];

outputMatrix = inputMatrix(index)';

end

A.2.6 inverseZigZag.m

function [ outputMatrix ] = inverseZigZag(inputMatrix)

%hard coded the inverse ZigZag

% 2. Write the MATLAB function idctmgr that implements the following

functionality:

% (a) It takes a matrix coeff of size 64 x Nblocks, where Nblocks is the

number of 8 x 8 blocks

% inside the image and reconstruct a luminance image.

% (b) For each given block b, the coefficients in the column coeff(:,b) are

used to reconstruct

% the block according to (3).

outputMatrix(1,1) = inputMatrix(1);

outputMatrix(1,2) = inputMatrix(2);

outputMatrix(2,1) = inputMatrix(3);

outputMatrix(3,1) = inputMatrix(4);

outputMatrix(2,2) = inputMatrix(5);

outputMatrix(1,3) = inputMatrix(6);

outputMatrix(1,4) = inputMatrix(7);

32

outputMatrix(2,3) = inputMatrix(8);

outputMatrix(3,2) = inputMatrix(9);

outputMatrix(4,1) = inputMatrix(10);

outputMatrix(5,1) = inputMatrix(11);

outputMatrix(4,2) = inputMatrix(12);

outputMatrix(3,3) = inputMatrix(13);

outputMatrix(2,4) = inputMatrix(14);

outputMatrix(1,5) = inputMatrix(15);

outputMatrix(1,6) = inputMatrix(16);

outputMatrix(2,5) = inputMatrix(17);

outputMatrix(3,4) = inputMatrix(18);

outputMatrix(4,3) = inputMatrix(19);

outputMatrix(5,2) = inputMatrix(20);

outputMatrix(6,1) = inputMatrix(21);

outputMatrix(7,1) = inputMatrix(22);

outputMatrix(6,2) = inputMatrix(23);

outputMatrix(5,3) = inputMatrix(24);

outputMatrix(4,4) = inputMatrix(25);

outputMatrix(3,5) = inputMatrix(26);

outputMatrix(2,6) = inputMatrix(27);

outputMatrix(1,7) = inputMatrix(28);

outputMatrix(1,8) = inputMatrix(29);

outputMatrix(2,7) = inputMatrix(30);

outputMatrix(3,6) = inputMatrix(31);

outputMatrix(4,5) = inputMatrix(32);

outputMatrix(5,4) = inputMatrix(33);

outputMatrix(6,3) = inputMatrix(34);

outputMatrix(7,2) = inputMatrix(35);

outputMatrix(8,1) = inputMatrix(36);

outputMatrix(8,2) = inputMatrix(37);

outputMatrix(7,3) = inputMatrix(38);

outputMatrix(6,4) = inputMatrix(39);

outputMatrix(5,5) = inputMatrix(40);

outputMatrix(4,6) = inputMatrix(41);

outputMatrix(3,7) = inputMatrix(42);

outputMatrix(2,8) = inputMatrix(43);

outputMatrix(3,8) = inputMatrix(44);

outputMatrix(4,7) = inputMatrix(45);

outputMatrix(5,6) = inputMatrix(46);

outputMatrix(6,5) = inputMatrix(47);

outputMatrix(7,4) = inputMatrix(48);

outputMatrix(8,3) = inputMatrix(49);

outputMatrix(8,4) = inputMatrix(50);

outputMatrix(7,5) = inputMatrix(51);

outputMatrix(6,6) = inputMatrix(52);

outputMatrix(5,7) = inputMatrix(53);

outputMatrix(4,8) = inputMatrix(54);

outputMatrix(5,8) = inputMatrix(55);

outputMatrix(6,7) = inputMatrix(56);

outputMatrix(7,6) = inputMatrix(57);

outputMatrix(8,5) = inputMatrix(58);

outputMatrix(8,6) = inputMatrix(59);

outputMatrix(7,7) = inputMatrix(60);

outputMatrix(6,8) = inputMatrix(61);

outputMatrix(7,8) = inputMatrix(62);

outputMatrix(8,7) = inputMatrix(63);

outputMatrix(8,8) = inputMatrix(64);

33

end

A.2.7 run_bits_value.m

function [symb] = run_bits_value(coefficientMatrix)

%implement run_bits_value

nZeros = zeros(64, length(coefficientMatrix));

nBits = zeros(64, length(coefficientMatrix));

zeros_n = 0;

for i = 1:length(coefficientMatrix)

for j = 1:64

value = coefficientMatrix(j,i);

nZeros(j,i) = -1;

if(value == 0)

zeros_n = zeros_n + 1;

end

if j == 1

nBits(j,i) = 12;

else

nBits(j,i) = 11;

end

symbMatrix = zeros(3*64, length(coefficientMatrix));

mZeros = zeros(64,length(coefficientMatrix));

mBits = zeros(64, length(coefficientMatrix));

mCoefficients = zeros(64,length(coefficientMatrix));

for i = 1:length(coefficientMatrix)

nPresent = 0;

for j = 1:64

if nZeros(j,i) ~= -1

nPresent = nPresent + 1;

mZeros(nPresent,i) = nZeros(j,i);

mBits(nPresent, i) = mBits(j,i);

mCoefficients(nPresent,i) = mCoefficients(j,i);

end

for i = 1:length(coefficientMatrix)

for j = 1:64

symbMatrix(1+3*(j-1),i) = mZeros(j,i);

symbMatrix(2+3*(j-1),i) = mBits(j,i);

symbMatrix(3+3*(j-1),i) = round(mCoefficients(j,i));

end

symbVector = zeros(3*64*length(coefficientMatrix),1);

numI = 1;

for i = 1:length(coefficientMatrix)

for j = 1:64*3 -1

if symbMatrix(j,i) == 0 && symbMatrix(j+1,i) == 0

symbVector(numI) = 0;

symbVector(numI + 1) = 0;

34

numI = numI +2;

break

else

symbVector(numI,1) = symbMatrix(j,i);

numI = numI + 1;

end

numI = 0;

for i = 1:length(symbVector) - 3

if symbVector(i) == 0 && symbVector(i+1) == 0 && symbVector(i+2) ==

0 && symbVector(i+3) == 0

break;

else

numI = numI +1;

end

symb = symbVector(1:numI);

end

A.2.8 inverse_run_bits_value.m

function [ coefficients ] = inverse_run_bits_value( symb )

%%(6)

blockN = 1;

start = zeros(4096,1);

start(1) = 1; %first block starts at 1

ends(4096) = length(symb); %ast block ends at the last symb element

for i = 3:length(symb)-1 %starting index of every block

if symb(i-2) == 0 && symb(i-1) == 0 && symb(i) == 0 %find EOB when three

consecutive zeros

start(blockN+1) = i;

blockN = blockN +1;

end

matComp = zeros(64,4096); %compressed matrix, doing decompression

for i = 1:4095

matComp(1:length(symb(start(i):start(i+1)-1)),i) =

symb(start(i):start(i+1)-1);

end

%handle last entry separately because of bound conditions

matComp(1:length(symb(start(4096):length(symb))),4096) =

symb(start(4096):length(symb));

matUncomp = zeros(64,4096); %uncompressed matrix

for i = 1:4096

ind = 0; %keep track of where the value is

for j = 1:3:61 %go in jumps of 3

zToSkip = matComp(j,i); %get the number of zeros that need to be

placed in between

nBits = matComp(j+1,i);

val = matComp(j+2,i); %get the value from the triplet

35

ind = ind+zToSkip+1; %new index to place val is the old

index+numZeros+1

matUncomp(ind,i) = val;

end

coefficients = matUncomp;